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Abstract 

We study the dynamics of one-dimensional solitons in the attractive and re- 
pulsive Bose-Einstein condensates (BECs) loaded into an optical lattice (OL), 
which is combined with an external parabolic potential. First, we demonstrate 
analytically that, in the repulsive BEC, where the soliton is of the gap type, 
its effective mass is negative. This gives rise to a prediction for the experi- 
ment: such a soliton cannot be not held by the usual parabolic trap, but it can 
be captured (performing harmonic oscillations) by an anti-trapping inverted 
parabolic potential. We also study the motion of the soliton a in long system, 
concluding that, in the cases of both the positive and negative mass, it moves 
freely, provided that its amplitude is below a certain critical value; above it, 
the soliton's velocity decreases due to the interaction with the OL. At a late 
stage, the damped motion becomes chaotic. We also investigate the evolution 
of a two-soliton pulse in the attractive model. The pulse generates a persistent 
breather, if its amplitude is not too large; otherwise, fusion into a single fun- 
damental soliton takes place. Collisions between two solitons captured in the 
parabolic trap or anti-trap are considered too. Depending on their amplitudes 
and phase difference, the solitons either perform stable oscillations, colliding 
indefinitely many times, or merge into a single soliton. Effects reported in this 
work for BECs can also be formulated for optical solitons in nonlinear photonic 
crystals. In particular, the capture of the negative-mass soliton in the anti-trap 
implies that a bright optical soliton in a self-defocusing medium with a periodic 
structure of the refractive index may be stable in an anti-waveguide. 
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1 Introduction 



An effective periodic potential in the form of an optical lattice (OL), created as 
an interference pattern between laser beams illuminating a Bose-Einstein con- 
densate (BEC), is a powerful tool that facilitates experimental studies of many 
dynamical properties of BECs. Among nonlinear dynamical perturbations that 
BECs may support are bright solitons, that were recently created in the exper- 
iment |2] • Direct observation of solitons in a BEC loaded in the OL has not yet 
been reported, although spitting of wave packets in a 87 Rb condensate, confined 
in a cylindrical trap with a weak longitudinal OL potential, which was observed 
in a very recent work [3] , can be explained as a manifestation of transient soliton 
dynamics. In anticipation of progress in the experiments, a lot of theoretical 
work has been done on the topic of BEC solitons in OLs ("solitons" are here 
realized as robust solitary waves, rather than as exact solutions of integrable 
equations). In particular, the use of OLs opens a way to observe gap solitons, 
that are expected to exist when the interaction between atoms in the BEC is 
repulsive (hence ordinary bright solitons are not possible) [5J E|. Moreover, 
while only one-dimensional (ID) solitons have been thus far observed in BECs, 
it was demonstrated that OLs may support multidimensional solitons of the gap 
type in the case of repulsion [5], and, which is plausibly more relevant to the 
experiment, they can also stabilize multidimensional 2D and even 3D solitons 
in an attractive BEC, which formally exist without the periodic potential, but 
are unstable because of the collapse [7j. Besides that, bright solitons with an 
internal vorticity can be stabilized too by the OL in the 2D case [7J|H]. OLs can 
also affect solitons in BECs in various other ways [§lllO|. 

Similar dynamical features are expected from optical solitons in nonlinear 
photonic crystals (NPCs) 111117115). In the latter case, a counterpart of the OL 
is periodic modulation of the local refractive index in the transverse direction(s) 
^1[7||H]. 2D solitons can be easily stabilized in the NPC JT] (in the uniform 
nonlinear medium, it would be the well-known Townes soliton, which is unstable 
due to the collapse [T3|L 

Thus far, most studies of solitons in BECs loaded in the OL did not system- 
atically consider their motion, nor two-soliton structures (very recently, mobility 
of discrete solitons, which represent the BEC solitons in OLs in the tight-binding 
approximation, was addressed, vis-a-vis the corresponding Peierls-Nabarro po- 
tential, in Ref. H). We aim to analyze these issues in the ID case. Several 
results will be reported. First, both ordinary solitons in the case of attraction, 
and gap solitons in the BEC with repulsion (the latter is more common in BECs 
|15|1 move across the OL freely if the soliton's amplitude (i.e., the number of 
atoms in it) does not exceed a certain critical value. Above that value, the 
soliton feels action of an effective braking force, induced by the Peierls-Nabarro 
potential, and eventually comes to a halt. It is relevant to mention that a 
similar property is known for moving lattice solitons in the discrete nonlinear 
Schrodinger (NLS) equation and related models: the solitons move freely unless 
they are too "heavy" ^JEl- F° r the solitons whose motion is hindered, we 
demonstrate dependence of the braking rate on the soliton's amplitude. At a 
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later stage of the braking, the soliton remains a coherent object, while its motion 
becomes chaotic. 

We also consider dynamics of a pulse which, in the absence of the OL po- 
tential, would give rise to a well-known two-soliton (breather solution) in the 
NLS equation (with attraction). We demonstrate that, if the initial amplitude 
of the two-soliton is below a critical value, it indeed develops into a persistent 
breather in presence of the OL. However, the two-soliton with the amplitude 
exceeding the critical value relaxes into a stationary (fundamental) soliton. 

It is known that the dispersion law for linear excitations, whose wavenumbcr 
is close to the edge of the Brillouin zone (BZ) induced by the OL potential, is 
characterized by a negative effective mass (curvature of the dispersion curve), 
unlike a vicinity of the BZ center, where the effective mass is positive. This 
feature was recently observed in an experiment JS], and it actually gives rise 
to the gap solitons in the case of the repulsive nonlinearity In this work, we 
will, first of all, demonstrate that the effective mass of the soliton proper may 
also be negative. A noteworthy feature of the analysis is that the underlying 
Gross-Pitaevskii (GP) equation, including the OL potential, does not conserve 
momentum, while equations providing for asymptotic description of the soliton 
dynamics conserve it. We provide for an explanation of this discrepancy. 

Then, we demonstrate that, having the negative mass, the gap soliton can- 
not be trapped by a usual external parabolic potential; however, it readily gets 
captured by an anti-trapping (inverted parabolic) potential (a possibility of con- 
fining a usual soliton in an inverted potential that rapidly oscillates in time was 
recently proposed in Ref. ^H], which, however, a completely different mecha- 
nism). This prediction suggests an experimental verification, especially in view 
of the fact that ordinary solitons in BECs were actually observed in inverted 
(expulsive) potentials 

Lastly, we also consider a state in which two identical solitons (with the pos- 
itive or negative mass) oscillate in opposite directions and periodically collide 
inside the, respectively, trap or anti-trap, applied on top of the OL. We demon- 
strate that there is a critical value of the amplitude of the solitons specific to this 
case, such that they merge into a single soliton (after one or several collisions) 
if the initial amplitude exceeds the critical value. Otherwise, the solitons keep 
to oscillate, colliding indefinitely many times. 

All the above-mentioned results can also be applied to optical solitons in 
NPCs. In particular, in the case of the self-defocusing nonlinearity (negative 
Kerr effect), the corresponding gap soliton will be stably confined in the NPC 
under the anti-trapping conditions, which, in optics, corresponds to an anti- 
waveguide I n fact , this suggests the first possibility to create a stable 
bright optical soliton in a nonlinear anti- waveguide. 

The rest of the paper is organized as follows. The model and analytical 
approximations for the solitons in it are presented in section 2. Numerical 
results for the motion of a soliton, and for two-solitons, are collected in section 
3. Trapping of the soliton by the parabolic potential is considered in section 4, 
and collisions between solitons oscillating in the trap or anti-trap are presented 
in section 5. Section 6 concludes the paper. 
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2 Formulation of the model and analytical ap- 
proximations for solitons 

2.1 Basic equations and momentum conservation 

The mean-field description of the BEC dynamics is based on the GP equation 
for the mean- field wave function ip in three dimensions |15|. 
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where m is the atomic mass, G = Airh 2 a/m, with a being the s-wave scattering 
length, and U is the external potential. As it was demonstrated in a number of 
works [20|, m the case of a quasi-lD ("cigar-shaped") configuration, Eq. (JTJ can 
be reduced to a ID equation. In the presence of the OL potential (the parabolic 
trap will be added later), the normalized equation for a ID mean-field wave 
function (f>(x,t) is [SHU: 

M _ ldV 

where a = +1 and — 1 for a > and a < 0, respectively (i.e., repulsion and 
attraction, respectively), the period of the OL is set to be 1, and — e is its 
strength. For the case of attraction, dynamics of solitons in this model was first 
investigated (in terms of nonlinear optics) in Ref. |12| . 

It is necessary to note that, while all the works published to date describe 
the gap solitons within the framework of the mean-field GP equation, effects of 
quantum depletion ( "evaporation" of atoms from the soliton wave function due 
to quantum fluctuations) may affect them similar to the way they were predicted 
to blur the notch of dark solitons in BECs as § a P solitons contain many 
notches [see. e.g., Fig. 1(b) below]. However, analysis of the quantum depletion 
is beyond the scope of the present work. 

First, we aim to give an explanation to the negativeness of the gap-soliton's 
mass in the case of u = +1 (repulsion). To this aim, we notice that the gap 
soliton may be represented by a combination of two terms, each being a product 
of a slowly varying amplitude u(x,t) or v(x,t) and rapidly varying carriers, 
exp (±iTrx): 

<f>(x, t) = (VE/2) [U{x, t)e™ x + V(x, t)e~™ x ] . (3) 

Substituting this into Eq. J2J i keeping only first derivatives of the slowly varying 
functions, and defining rescaled variables 

r ee {e/2)t, £ = {e/2tt)x (4) 

lead to the standard model (J5J that gives rise to gap solitons: 

.8U . ,dU /1,„ |2 lT , l2 



/ — = --— — - \a\4>\ 2 - ecos(27rx)] 0, (2) 



.dV ,8V fl 



i— 1- 



8t d£ \2 



\V\ 2 + \U\ 2 ) V + U = 0. (5) 
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Note that both the underlying GP equation (J2J and the asymptotic gap- 
soliton model (JSJ conserve the norm (number of atoms). In the former case, it 
is 

/+oo 
<j) 2 (x)dx. (6) 
-oo 

The substitution of the waveform (j3J into this expression and neglecting terms 
like J_°° U(x)V*(x)e 2mx dx (they are exponentially small, being integrals of 
products of rapidly oscillating functions and slowly varying ones) yield N = 
(tt/2) [|/7(0| 2 + \V{0\ 2 ] dx, which is indeed the norm conserved by Eqs. 

Q 

On the other hand, Eqs. JSJ conserve the momentum, 

f +oc fdU* dV* \ 



(the asterisk stands for the complex conjugation), while the momentum of Eq. 

P GP = 2ij -^-cj> dx, (8) 

is not conserved in the presence of the periodic potential, but rather evolves 
according to the equation 

dP G - r+cc 



dt 



r+oo 

:P = -Aire / |</)(j:)| 2 sin(27ra;)dx, (9) 

J — oo 



(integration by parts was used to simplify the expression). Making use of Eqs. 
121 and |0J and dropping exponentially small terms of the above-mentioned 
type, one can cast Eq. © in the form 

= 2n 2 i [ + °° (U*V - UV*) d£. (10) 

To explain the apparent contradiction between the conservation of P and 
non-conservation of Pqp, we note that the substitution of Eqs. 10 and (0J into 
Eq. (JSJ) yields, aside from exponentially small terms, an expression 



Pgp = 



d^. 



Inserting this, and a relation which is a consequence of Eqs. (JSJ), 

J I- + CO r + oo 

"T / (\U\ 2 -\V\ 2 )d£ = 2i / (U*V-UV*)dt (11) 

^ J —oo J — oo 

into Eq. © results in dP/dr = 0. 
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2.2 The mass of gap solitons 



An analytical solution to Eqs. (JSJ corresponding to a gap soliton moving at a 
velocity v, which belongs to the interval — 1 < u < +1, is well known |23- In 
the case of a = —1 (the attractive BEC), it is 



U = 



V 



1 



2(1-1/ 



>V 4 
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x exp 



1-1! 



■ V 
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1 ) tan 1 ( tan — ■ tanh (3 sin ( 



iT cos 6 



9\l/4 . 

v sin ( 



3 - v 2 ^cosh (25 sin 9) + cos 



x exp 



4d 



1 tan 



tan — • tanh (3 sin ( 



iT cos 6 



(12) 



— 1/2 

v 2 ) (r — v£), and 8, which takes 



where S = (l - v 2 ) ' (£-vt). t =(1 
values < 6* < tt, determines the amplitude and width of the soliton. For the 
soliton l|12l) . the momentum Q can be calculated in an exact form: 



3- 



7- 



3 — 1!" 



(sin 8 - 0cos0) + 8 cos 6» 



(13) 



The expression (|13|l simplifies for slow solitons (v 2 <C 1), making it possible to 
define the mass of the slow gap soliton as the momentum divided by the velocity, 



m (GS) = P 

'"■sol — r » 



,i/u= (8/9) (7sin6»-46»cos6') 



(14) 



Note that this mass is positive for all values of 8. 

In the case of the repulsive BEC, a = +1, we obtain Eqs. © with the 
opposite sign in front of the cubic terms. To cast the equations into the standard 
form to which the soliton solution (|12f> pertains, we apply complex conjugation 
to them, and define U = U*,V = — V* . Then, the gap-soliton's momentum 
P defined by Eq. Q with U and V replaced by U and V is related to the 
velocity exactly the same way as above, i.e., as per Eq. I|13fl . However, the 
proper momentum P is still defined by Eq. 10 in terms of the original fields 
U and V, hence P = —P. Thus, the effective mass in the repulsive model is 
precisely -m[j S ' , being always negative. 



2.3 Description in terms of the Bloch functions 

To continue the analysis, we now return to the underlying GP equation PJl. 
Another approach to solitons in the weakly nonlinear case starts with the linear 
Bloch functions that obey the Mathieu equation, 

EF(x) = (1/2) F" + e cos (2ttx) F. (15) 
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Here E is the chemical potential, and the solution is quasi-periodic, i.e., F(x) — 
e lkx f(x), where A; is a wavenumber, and the function f(x) is periodic, f(x + l) = 
f(x). The solution determines the corresponding band structure, E = E(k), 
the points k = and k — tt being, respectively, the center and edge of the first 
Brillouin zone (BZ). 

Then, an approximate solution to the underlying GP equation (J2Jl for a 
small-amplitude broad soliton is looked for as 

<p{x, t) = e- lEt F{x)<S>{x, i), (16) 

with a slowly varying envelope function $(a;), cf. Eq. @. As the norm of 
the solution (the number of atoms in the BEC), N = J_ \4>\ 2 dx, approaches 
zero, the ordinary soliton solution in the attractive model and its gap-soliton 
counterpart in the repulsive one approach, respectively, the linear Bloch function 
with k = and k = tt. 

In the general case, the slowly varying amplitude obeys an asymptotic NLS 
equation, that can be derived by substituting the ansatz (|TTfl) into Eq. @ 5 : 

.d$ 1 _i9 2 $ ,_. a _ 

t at = -2 m ^ + tra ^ 9 - (17) 

Here, the effective mass m c g for linear excitations is determined from the cur- 
vature of the band structure of the linear Bloch states as m~ s = E"(k), and 
g = J Q \F(x)\ 4: dx/ J Q \F(x)\ 2 dx (recall 1 is the normalized period of the OL po- 
tential). Notice that Eq. Ijl7|l conserves the momentum, unlike the underlying 
GP equation J5J. This difference can be explained as it was done above, see 
Eqs. - GBJ. 

In the case of m e ff<7 < 0, Eq. (|17|) has soliton solutions. Their exact form is 
$(x) = A exp [i (kx - AE ■ t)] sech \Ay/g |m e g | (x - vt)j (18) 
with the wavenumber k, velocity 

v = k/m cS , (19) 

an arbitrary amplitude A, and 

AE = -i 5 A 2 sgn (m cff ) + ^m cB v 2 . (20) 

Close to the BZ center, i.e., for small k, the Bloch function may be reasonably 
approximated by a combination of three relevant harmonics, 

F(x) = c\ exp(ikx) + ci exp (i(k + 2n)x) + C3 exp (i(fc — 2tt)x) . (21) 

Substitution of this into Eq. I|15(l yields, in the first approximation, E(k) = 



^(center) + fc2/ ^^center)^ ^ ^ 



^(center) = ^ _ ~ ^ ^ 
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and an expression for the effective mass which is indeed positive, 



in 



(center) 



27r 4 + e 2 +7r 2 v / 47r 4 + 2e 2 



fOvr 4 + e 2 - 37r 2 V4vr 4 + 2e 2 



(23) 



Other coefficients in Eqs. Ij2f(l and (|20|) are found to be c\ — f , C2 = C3 = —Eq/s, 
and 5 = (1 + 12c 2 , + 6c 4 )/(l + 2cf). Then, with regard to Eqs. JTHJ and flty . 
the full lowest-order approximation (corresponding to k — 0) for the localized 
state in the attractive model is 



4>{x,t) = A exp 



(center) 

im^j vx 



E 



(center) 



AE)t 



[1 + 2c2 cos(27rx)] sech A 




(x -vt)\, (24) 



where AE is given by Eq. I|20(l . We stress that this soliton exists only in the 
attractive model, with a = — 1. 

In a similar way, one can work out an approximation for the linear Bloch 
function at the BZ edge, i.e., near k = n, as a combination of two harmonics: 



F(x) = c\ exp(ikx) + C2 exp (i(k — 2t:)x) , 

cf. Eq. ©. The substitution of this into Eq. JUJ yields E(k) = E ( 
(k~Tr) 2 /(2m^ ff dgo) ), where 



(25) 



(edge) 



E, 



(edge) 



(^ 2 ±N)A 



(edge) 



(26) 

'" cff " leliL 2 ' (27) 

and g = 3/4. The signs + and — corresponds to the second and first Bloch 
bands, respectively. As is seen, the effective mass (|27|l for the linear excitations 
near the BZ edge for the first Bloch band is negative if |e| < 2-7T 2 (below, we 
present typical numerical results for e = 5, which satisfies this condition), and 
the effective mass for the second Bloch band is positive. These two possibilities 
correlate to the conclusion formulated above within the framework of the ap- 
proximation based on Eqs. © and JSJ, that there may exist gap solitons with 
positive and negative masses [that approximation definitely corresponds to the 
case of the weak OL, i.e., |e| < 27r 2 , in terms of Eq. l(77)l]. 

The full form of lowest-order approximation (corresponding to k — ir) for 
the soliton carried by the Bloch wave near the BZ edge in the first Bloch band 
(the one with the negative effective mass) follows from Eq. l(T%|) [cf . Eq. (|2"4fl ] : 



4>(x) = A exp 



(edge) 

jrn^j vx 



e) 



x cos(7rx) sech A\ g 



(edge) 



AE) t 



(x — vt) 



(28) 



As it follows from Eq. I|17l) . this soliton exists in the repulsive model (er = +1) if 
|e| < 2-7T 2 . On the other hand, the soliton carried by the Bloch wave near the BZ 
edge in the second band exists in the attractive model, since it has rn^ 80 ' > 0. 
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3 Standing and moving solitons 



3.1 Zero- velocity solitons 

Stationary soliton solutions can be easily found in the form o{<f>(x, t) — 4>{x) exp (—iEt), 
applying the numerical shooting method to the stationary version of the un- 
derlying GP equation @. Figures 1(a) and 1(b) display, respectively, typical 
examples of the localized solution for the attractive (a = —1) and repulsive 
(cr = +1) model with the OL potential (e = 5). In the former case, the analyti- 
cal approximation (1221) with e = 5 yields the value Eq = —1.194 of the chemical 
potential at the BZ center, and in the latter case, Eq. I|2t)|) with e — 5 yields 
E = -7. 565 at the BZ edge. 



0(x) (a) 0(x) (b) 




-10 -5 5 10 -10 -5 5 10 

X X 



Figure 1: Typical examples of numerically found and analytically predicted 
(solid and dashed curves, respectively) stable zero- velocity solitons in Eq. J5J 
with attraction, a = — 1 (a), and repulsione, a — +1 (b). The amplitude of the 
optical-lattice potential is e = 5. 

In the example shown in Fig. 1(a), the soliton's amplitude is A = 0.5 and 
its norm is N = 0.354. In Fig. 1(b), the amplitude is A = 1, and the norm 
is N = 0.873. In both parts of Fig. 1, the dashed profiles show the analytical 
predictions 1241) and (|28ll , respectively, corresponding to the same values of the 
amplitude. The comparison with the analytical approximations makes it obvious 
that they are very accurate in both cases. At other values of the parameters, 
the agreement between the numerically found and analytically predicted shapes 
of the solitons is, generally, as good as in these examples. 

It is necessary to stress that, alongside broad solitons extending to many 
periods of the OL, a generic example of which is shown in Fig. 1(a), stable 
narrow solitons, which are confined, essentially, to a single OL cell, also exist 
in the attractive model. They were found in Ref. (T2|. Quite similarly, the 2D 
and 3D counterparts of the ID equation (J2J with attraction also support both 
narrow ( "single-cell" ) and broad ( "multi-cell" ) stable solitons . 
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Figure 2: Results of numerical simulations of the moving soliton generated by 
the initial condition (|29|l with k = 0.2. The panels (a) and (b) here, as well as 
in Figs. 4 and 5 below, and Figs. 6 and 8 show the time evolution of |</>(x)|. In 
the panel (a), the initial amplitude of the soliton is A = 0.5, while panels (b) 
and (c) pertain to A = 0.75. The panel (c) displays a relatively early stage of 
the motion of the soliton's center in the case of braking, the dashed curve being 
the fit to Eq. 1(30(1 ■ The panel (d) shows the braking rate, defined as in Eq. 
1)30(1 . vs. the soliton's amplitude. 



3.2 Moving solitons in the attractive model 

The next step is to simulate moving solitons that were also predicted above. 
Starting with the attractive model, in Fig. 2 we display a numerically found 
moving soliton, which is generated, in the attractive model, by the analytical 
waveform 1(24(1 with e = 5, 

<t>o(x) = A e lk(x - L/2) [0.877 + 0.216 cos (2w(x - L/2))] sech (A(x - L/2)) , 

(29) 

taken as the initial condition for Eq. J2J. In this case, periodic boundary con- 
dition were used with the spatial period L = 50 (results of simulations with the 
parabolic trapping potential will be presented below). The initial wavenumbcr 
in Eq. 1(29(1 is k — 0.2; the amplitude is A = 0.5 in the case shown in Fig. 2(a), 
and A = 0.75 in Fig. 2(b) (the value of k — 0.2 is formally incompatible with 
the system's period, L = 50, but this mismatch plays no practical role, as the 
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soliton's size, l so \ ~ A^ 1 , inside which the phase field carrying the wavenumbcr 
was actually created, is much smaller than L). 

For A = 0.5, a stable soliton steadily moving at the velocity v = 0.175 was 
found in the simulations, see Fig. 2(a) (the velocity was identified as that of 
the soliton's peak). Following Eq. (|19fl . the effective mass corresponding to it 
was calculated as k/v, with the wavenumber k of the initial configuration (|29l) . 
The thus found effective mass 1.143 is very close to the analytical value 1.13 
produced by Eq. J23) . 

For A = 0.75, the simulations produce an altogether different result: the 
soliton's velocity slowly decreases with time, see Figs. 2(b) and 2(c). At a 
relatively early stage of the evolution, the time dependence of the velocity may 
be well fitted by 

v(t) =v exp(-jt), with v = 0.17,7 = 0.0025. (30) 

To illustrate this, the law of motion X p {t) = 25+0.17Q1 - exp(-0.0025t)] /0.0025, 
which is the temporal integration of Eq. i|30[l . is depicted in Fig. 2(c) by a dashed 
curve, and the actual trajectory of the soliton's center is shown by a continuous 
line. 

Figure 2(d) displays the braking rate 7, which is defined as in Eq. l|37)|) . i.e., 
by fitting the initial law of motion to the exponential decay of the velocity, as 
a function of the soliton's amplitude A. This dependence was obtained from 
systematic simulations of Eq. with the fixed OL strength, e = 5, and 
different initial conditions. The braking sets in at the critical value of the 
amplitude, A CI — 0.72. It is relevant to mention that a qualitatively similar 
effect was observed, in various forms, in simulations of the motion of a discrete 
soliton in the NLS lattice equation: if the amplitude of the moving soliton 
exceeds a critical value, it quickly gets trapped by the lattice |17| . otherwise it 
moves freely [To) . 




t 



Figure 3: The decay of the soliton's momentum, in the same case as shown in 
Figs. 2(b) and 2(c), on a longer time scale. The dashed curve is the fit to Eq. 
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At a later stage of the evolution, the soliton's motion law strongly deviates 
from the exponential relaxation assumed in Eq. (|30() . To illustrate this, in Fig. 3 
we display the evolution of the soliton's momentum per atom, which is computed 
from the numerical data as p = i J^°° {dip* /dx) <f>(x)dx/ \<f>(x)\ 2 dx [cf. Eq. 
iJSJ], at a moderately long time scale. The initial exponential decay, p(t) « 
0.168 exp(— 0.0028i), with the decay rate very close to that in Eq. Ipffi)l. changes 
into an erratic motion at a later stage, with the deceleration sometimes changed 
by acceleration. This picture suggests that dynamical chaos may develop in the 
motion of the soliton. To further investigate the issue, in Fig. 4 we display 
continuation of the motion to extremely long time. The figure, and, especially, 
the continuous power spectrum in panel (c) clearly demonstrate that chaotic 
motion sets in indeed. However, the soliton's shape remains strictly coherent at 
all times, hence that the mean-field approximation based on the GP equation 
remains relevant in this case, unlike situations when the wave field as a whole 
becomes chaotic, and the BEC disintegrates, making the GP equation irrelevant 

E2- 




Figure 4: The continuation of the Fig. 3 to an extremely long time scale. The 
panels (a) and (b) display the time dependences of the coordinate of the soliton's 
center and momentum per atom, and the panel (c) is the power spectrum of the 
latter in the double-logarithmic scale. 



3.3 Relaxation of a two-soliton in the attractive model 

It is well known that the usual self- focusing NLS equation, i.e., Eq. @ with 
g = —1 and e = 0, gives rise not only to fundamental solitons, but also to their 
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higher-order counterparts, the simplest one being the two-soliton (a breather), 
which is generated by the initial condition with the double amplitude. We sim- 
ulated evolution of similar configurations in the presence of the OL potential, 
taking, for example, the initial condition as in Eq. I|29[) but with the amplitude 
2A instead of A [while sech(A(x — L/2)) was not altered]. If the amplitude is 
smaller than some critical value, this configuration indeed evolves into a persis- 
tent breather, see an example in Fig. 5(a) for A = 0.16. However, a "heavier" 
initial two-soliton configuration relaxes into a stationary fundamental soliton, 
like in Fig. 5(b) in the case of A = 0.3. Figure 5(c) additionally displays the 
evolution of the pulse's amplitude in the latter case. 




Figure 5: Evolution of the two-soliton, initiated by the initial configuration 1|29|1 
with the double amplitude: (a) A — 0.16; (b) and (c): A — 0.3. A persistent 
breather appears in the former case, while in the latter case the pulse relaxes 
(without conspicuous radiation loss) into a fundamental soliton. 

It is relevant to mention that the attractive model with the OL potential 
may have stable static two-soliton (or multi-soliton) solutions of a different type, 
when two narrow solitons are trapped in two adjacent wells of the OL potential. 
Such states and their stability limits were found in Ref. ^2] in terms of a 
nonlinear-optical model. 

3.4 Moving solitons in the repulsive model 

Motion of solitons in the repulsive GP equation J2Jl, with a = +1 , was sys- 
tematically investigated too. Again fixing e = 5, we take, as a generic example, 
the initial condition in the form suggested by the analytical approximation (|28() , 
i.e., 

(j> {x) = 1.985 A e lk{x - L/2) cos (tt(x - L/2)) sech (A(x - L/2)) , (31) 
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cf. Eq. (|29l) . The outcome of the simulation, corresponding to k = 0.1 in Eq. 

is shown in Fig. 6, which displays the evolution of \4>(%)\ for two different 
amplitudes, (a) A = 0.30 and (b) A = 0.55. In the former case, the gap soliton 
moves steadily at the velocity v — —0.274. The effective mass corresponding 
to this soliton, as determined from the numerical data according to the same 
definition as above, m e s = k/v, is —0.365, which is close enough to the analytical 
value -0.34 produced by Eq. (27)l. 



(a) 




0.48 0.49 0.5 0.51 52 53 54 
A 



Figure 6: The panels (a), (b) and (c) show the same as (a), (b) and (d) in Fig. 
2, but for the case when the soliton is launched in the repulsive model, using 
the initial condition H31(l with k = 0.1. The panels (a) and (b) correspond to 
the amplitudes A — 0.30 and 0.55, respectively. 

Figure 6(b) shows that a heavier soliton, with A = 0.55, does not move 
steadily. Instead, it starts to brake, like heavy solitons in the attractive model, 
see Fig. 2. The velocity decays exponentially at the initial stage of the evolution, 
like in Eq. J3UJ), the corresponding damping rate 7 being displayed in Fig. 6(c) 
versus A. As well as in the attractive model, the transition from the steady 
motion to the braking occurs at a critical value of the soliton's amplitude, which 
is A cr ss 0.48 in this case. 



4 Solitons in the parabolic trap or anti-trap 

Any experimental setup in BEC uses a parabolic trap \T5\ (bright ID solitons 
were observed in an expulsive potential created by an inverted trap 2 ). This 
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means that the underlying GP equation should be extended to the form 



d<j> 



ld 2 <p 



cr\4>\ - e cos (2irx) + B I x 



where B is the strength of the trap, and its center is set at x = L/2. 



(32) 




Figure 7: Motion of the soliton in the weak parabolic potential: (a) the trap 
with B = 0.0001; (b) the inverted potential (anti-trap) with B = -0.0001. The 
strength of the optical lattice is e = 5. In the cases shown in (a) and (b), the 
solitons were set in motion by lending them the initial velocities vq = 0.172 and 
—0.270, respectively. 

Figure 7(a) shows the evolution of \<j>(x)\ in the attractive model (a — —1) 
with e = 5 and a weak trap, B = 0.0001, as obtained from numerical simulations 
of Eq. IML't . The initial condition is again taken in the form of Eq. I|29|) , which 
is suggested by the analytical approximation l|24(l . with A = 0.5 and k — 0.2. In 
this case, the soliton performs harmonic oscillations in the trap, and the motion 
of the soliton's center-of-mass coordinate, which we define as 

/OO POC 
x\<f>\ 2 dx/ / \$\ 2 dx, (33) 
-OO J —OG 

is well approximated by the expression 

X p (t) = 25 + 13sin(0.0132 t). (34) 

To proceed with an analytical approach to the soliton's motion in the parabolic 
potential, we insert the same ansatz Qltty . which was used above, into Eq. (|32|) . 
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and derive a modified version of Eq. i|17p. 

i<$> t = - (1/2) m cff f H + [ag\$\ 2 + B(x - L/2) 2 ] $. (35) 
Then, the following relation is an exact corollary of Eq. (|35[1 . 

m cS X p = -2B(X p - L/2), (36) 

where the overdot stands for the time derivative. In the model with attraction, 
where the effective mass is positive, Eq. (|36|l predicts harmonic oscillations of 
the soliton with the frequency uj\ 10 = ^/2B/m c s . For instance, in the example 
shown in Fig. 7(a), Eq. I|23|l with e = 5 yields m c s — 1.13, hence, with 
B = 0.0001, we predict the frequency uJho — 0.0133, which is to be compared 
with the numerically found value lo num = 0.0132 in Eq. (|34|l . 

Figure 7(b) shows the evolution of |</>(a;)| in the repulsive model with the 
inverted parabolic potential, B = —0.0001. The initial condition is again taken 
in the form of Eq. H31(l . with A = 0.3 and k = 0.1. It is clearly observed that 
the negative-mass soliton is trapped by the anti-trapping potential. The motion 
of the soliton's center is fitted by 

X p {t) = 25 - 11.5sin(0.0235 t), (37) 

cf. Eq. H34fl . In this case, the analytical approximation H27fl with e = 5 yields 
m e ff = —0.339, hence the corresponding harmonic-oscillation frequency with 
B = —0.0001 is uj\ w = ^/2B/m c g = 0.0243, which is close to the numerically 
found value w nU m that appears in Eq. (|3*7j) . 

We have checked too that, in the case of the usual trap (B > 0), the negative- 
mass soliton in the repulsive model is indeed expelled from the trap region, as it 
should be expected. In that case, its observed motion again complies well with 
Eq. p)|). 

Lastly, we have also found that, as well as in the absence of the parabolic 
potential (see above), the motion of the soliton is subject to braking if its am- 
plitude is too large. In that case, the positive- and negative-mass solitons per- 
form damped oscillations in the trapping or anti-trapping potential, respectively. 
Typical examples are shown in Fig. 8, which displays the motion of the soliton's 
center as found from the simulations initiated by the same initial conditions l|29|l 
and H31(l as above, but with larger amplitudes: A = 0.75 in (a), and A = 0.52 
in (b). In these two cases, the damped oscillations are well approximated by 
analytical expressions 

X p (t) = 25 + 12.5 sin(0.0132i) exp(-0.0010 t) (38) 

[the dashed curve in Fig. 8 (a)], and 

X p [t) = 25-10 sin(0.0202 t) exp(-0.0006 t) (39) 

[the dashed curve in Fig. 8 (b)]. 

At other values of the parameters, the motion of the positive- and negative- 
mass solitons in the trapping and anti-trapping potentials is the same as shown 
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Figure 8: Numerically found trajectories of the soliton's center in the presence 
of the harmonic trap or anti-trap, in the case when the soliton's amplitude 
exceeds the braking threshold. The GP equation (|32|l was simulated, with 
a = —1 and B = 0.001 (attraction and normal trap) in (a), and with a = +1 
and B = —0.001 (repulsion and anti-trap) in (b). The dashed curves show the 
respective analytical fits (|38|l and i|39|) . 



in the above examples. In all the cases when the braking did not set in, the 
numerically simulated motion (both oscillations and expulsion, depending on the 
signs of the effective mass and the potential) was found to be in good agreement 
with the predictions based on Eq. (|3fcil> - 

The trapping of the gap soliton by the inverted parabolic potential in the 
repulsive BEC loaded into an OL is a challenge for experimental verification. We 
stress that, while the gap soliton should be trapped by the inverted potential, 
this is definitely not expected for a plain broad distribution of the atomic density, 
in which case the anti-trapping potential will try to expel atoms in the usual 
way (of course, this trend may be countered by trapping the atoms in the OL). 
Thus, the trapping in the inverted parabolic potential is a specifically nonlinear 
dynamical effect, endemic to the negative-mass solitons. 



5 Collisions between oscillating solitons 

If the direct or inverted parabolic potential retains the solitons, it is quite easy 
to create an initial configuration with two well-separated identical solitons, and 
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lend them opposite initial velocities (or just zero velocities, see below). Then, 
they will start to oscillate in opposite directions, and will thus have to collide 
periodically. This setting is quite convenient to test the character of collisions 
between solitons. 




Figure 9: Collisions between positive-mass (a) and negative-mass (b,c) solitons 
which oscillate, respectively, in the trapping or anti-trapping parabolic potential 
with the strength B = ±0.0001. The difference between the cases shown in (b) 
and (c) is in the initial phase difference Atpo between the negative-mass solitons: 
Aifo = t in (b), and Atpo = in (c). 

We simulated this situation in the attractive model, creating an initial config- 
uration (/)q (x) which was a linear superposition of the analytical approximations 
(|24[) for two solitons. As a typical example, in Fig. 9(a) we display the case 
when centers of the two solitons with equal amplitudes A = 0.3 are set at the 
points x\ = L/2 — 13 and x 2 — L/2 + 13, so that 

<f> (x) = A{[0.877 + 0.216cos27r(a;-xi)]sech(A(x-xi)) 

+ [0.877 + 0.216 cos 2n(x - x 2 )]sech{A(x - x 2 ))} ■ (40) 

Note that the initial wavenumbers in this expression are set equal to zero [cf. 
Eq. I|29(l ]. hence the initial velocities of the two solitons are zero too; however, 
being in-phase, the solitons attract each other, which initiates their oscillatory 
motion. Figure 9(a) clearly demonstrates (which is also corroborated by ex- 
tending the simulations to much longer times) that the positive-mass solitons 
collide repeatedly, each time passing through each other, and remain completely 
stable. 

A typical example of a gap-soliton pair that survives indefinitely many colli- 
sions in the repulsive BEC placed in the inverted parabolic potential is displayed 
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in Fig. 9(b). The two solitons were initially taken as a linear superposition of 
two analytical waveforms Q28[l. again with equal amplitudes A = 0.3, centers 
placed at the points x\ = L/2 — 13 and x 2 = L/2 + 13,and zero initial wavenum- 
bers: 

4>o(x) = 1.985A [cos (ir(x — x±)) sech (A(x — xi)) — cos (ir(x — x^)) sech (A(x — x 2 

(41) 

Unlike the initial condition (|29|l which was used in the attractive model, the 
initial phase difference between the solitons in the configuration (|40|l is Atpo = 
hence they repel each other. In this case too, the solitons perform indefinitely 
many oscillations, colliding repeatedly without developing any instability. Note 
that, in contrast with the case shown in Fig. 9(a), this time the colliding solitons 
do not pass through each other, but rather bounce, which is explained by the 
fact that they keep the phase difference 7r. 

However, the solitons do not always collide elastically. Another generic out- 
come is fusion of two in-phase solitons into a single one. For example, in the 
repulsive model with the anti-trapping parabolic potential and the same initial 
configuration as in Eq. Ij41|l . but this time with Aipo = 0, 

4>o(x) — 1.985^4 [cos (tt(x — xi)) sech (A(x — xi)) + cos (tt(x — x 2 )) sech (A(x — x 2 

we have found that the first collision leads to merger of the solitons into one 
pulse, which is accompanied by some radiation loss, see Fig. 9(c). The resul- 
tant single soliton then performs harmonic oscillations with a small amplitude. 
Qualitatively, the fusion may be explained the following way: in the course of 
the collision, overlapping between the in-phase solitons leads to increase of the 
amplitude past the threshold at which the braking sets in (see above). 

Further simulations demonstrate that, in the case shown in Figs. 9(b) and 
9(c), the colliding solitons eventually merge into one if the initial phase difference 
between them belongs to the interval |A<^o| < 0.87T. It is noteworthy that, when 
Aipo is close to the critical value (0.87T, in this example), the fusion happens 
after several collisions. 

In the attractive model, the merger induced by the collisions was observed 
too, if the solitons' amplitude was large enough. On the other hand, in both 
models (repulsive and attractive) merger of solitons with Aipo = ir has not been 
observed. We stress that (as well as in the case of the discrete NLS equation 
|17|1 solitons with very large amplitudes cannot collide at all, as they will be 
quickly halted by the braking force induced by the OL. 

6 Conclusion 

In this work, we aimed to investigate, in a systematic way, motion of ID solitons 
in the attractive and repulsive Bose-Einstein condensates (BECs) loaded into an 
optical lattice (OL) and subject to the action of an external parabolic potential. 
First of all, we have demonstrated, in two different analytical forms, that, in the 
repulsive BEC, when the soliton is of the gap type, its effective mass is negative. 
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This gives rise to the prediction that calls for experimental verification: in 
the latter case, the soliton cannot be confined by the usual trapping parabolic 
potential, but it can be held by the anti-trapping inverted potential. We have 
also investigated the motion of the solitons in long systems, concluding that, 
in both cases of the positive and negative mass, the soliton moves unhindered 
if its amplitude is below a certain threshold; above the threshold, the soliton 
is braked by its interaction with the OL. At a late stage of the braking, the 
damped motion becomes chaotic. 

We have also investigated the evolution of the two-soliton state in the attrac- 
tive model. It was found that the two-soliton generates a persistent breather, 
if its amplitude is not too large; otherwise, the two-soliton suffers fusion into a 
single fundamental soliton. 

Lastly, we have considered collisions of two solitons trapped in the parabolic 
potential. Depending on their amplitudes and phase difference Aip, the solitons 
perform stable oscillations, colliding indefinitely many times, or they merge into 
a single soliton, after one or several collisions. However, the merger does not 
occur in the case of Aip = tt. 

As concerns values of physical parameters at which the phenomena predicted 
in this work can be observed in the experiment, they are not different from those 
for which other soliton effects in OLs were recently predicted [SI E| (the number 
of atoms in the soliton may be in the range of 10 3 — 10 4 , the number of the 
OL periods inside the parabolic trap or anti-trap can easily be ~ 100 or larger, 
etc.). Thus, the observation of the gap-soliton's confinement in the anti-trapping 
potential and other dynamical phenomena predicted in this work seems quite 
feasible. Besides that, all the predictions may also be translated into those for 
solitons in nonlinear photonic crystals. In particular, the capture of the gap 
soliton by the inverted potential implies that a stable bright optical soliton is 
possible in a self-defocusing nonlinear anti- waveguide, with the refractive index 
periodically modulated in the transverse direction. 

Results reported in this paper can be extended to the two-dimensional (2D) 
case. In particular, we have demonstrated that a 2D gap soliton in a repulsive 
BEC loaded in a 2D optical lattice may be trapped by an inverted isotropic 
harmonic potential. In that case, the trapped soliton demonstrates various 
modes of stable motion, including oscillations with periodic passage through the 
center, and circular motion around the center. These results will be reported in 
detail elsewhere. 
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